Estimating $\pi$ by Sampling Points

By Evgenia "Jenny" Nitishinskaya and Delaney Granizo-Mackenzie

Notebook released under the Creative Commons Attribution 4.0 License.


A stochastic way to estimate the value of $\pi$ is to sample points from a square area. Some of the points will fall within the area of a circle as defined by $x^2 + y^2 = 1$, we count what percentage all points fall within this area, which allows us to estimate the area of the circle and therefore $\pi$.


In [34]:
# Import libraries
import math
import numpy as np
import matplotlib.pyplot as plt

In [45]:
in_circle = 0
outside_circle = 0

n = 10 ** 4

# Draw many random points
X = np.random.rand(n)
Y = np.random.rand(n)

for i in range(n):
    
    if X[i]**2 + Y[i]**2 > 1:
        outside_circle += 1
    else:
        in_circle += 1

area_of_quarter_circle = float(in_circle)/(in_circle + outside_circle)
pi_estimate = area_of_circle = area_of_quarter_circle * 4

In [46]:
pi_estimate


Out[46]:
3.1328

We can visualize the process to see how it works.


In [47]:
# Plot a circle for reference
circle1=plt.Circle((0,0),1,color='r', fill=False, lw=2)
fig = plt.gcf()
fig.gca().add_artist(circle1)
# Set the axis limits so the circle doesn't look skewed
plt.xlim((0, 1.8))
plt.ylim((0, 1.2))
plt.scatter(X, Y)


Out[47]:
<matplotlib.collections.PathCollection at 0x7f331b408d10>

Finally, let's see how our estimate gets better as we increase $n$. We'll do this by computing the estimate for $\pi$ at each step and plotting that estimate to see how it converges.


In [39]:
in_circle = 0
outside_circle = 0

n = 10 ** 3

# Draw many random points
X = np.random.rand(n)
Y = np.random.rand(n)

# Make a new array
pi = np.ndarray(n)

for i in range(n):
    
    if X[i]**2 + Y[i]**2 > 1:
        outside_circle += 1
    else:
        in_circle += 1

    area_of_quarter_circle = float(in_circle)/(in_circle + outside_circle)
    pi_estimate = area_of_circle = area_of_quarter_circle * 4
    
    pi[i] = pi_estimate
    
plt.plot(range(n), pi)
plt.xlabel('n')
plt.ylabel('pi estimate')
plt.plot(range(n), [math.pi] * n)


Out[39]:
[<matplotlib.lines.Line2D at 0x7f331b80c6d0>]